Efficient generation of realistic model systems of amorphous silica 
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We used classical molecular dynamics and the van Beest Kramer van Santen (BKS) potential 
to generate small model systems of amorphous silica. We further optimized the classically equili- 
brated configurations using plane wave based density functional theory and a generalized gradient 
(GGA) approximation. Within ab initio treatment we showed that both geometry optimization 
and Car-Parrinello annealing lead to the same final configurations but the CPU time required for 
the geometry optimization to reach convergence is one fifth of the time needed by a Car-Parrinello 
annealing. In addition during the optimization or the annealing no substantial change occurs on the 
topology acquired by the vitreous silica at the end of the classical quenching protocol. Structural 
and electronic properties have been calculated and compared to experiments. 
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I. INTRODUCTION 

Disordered forms of silicon dioxide {Si02) glasses play 
a key role in electronic device applications like semicon- 
ductor devices and optical fibers. Amorphous silica has 
a structure consisting of a continuous random network 
of corner-sharing 5*^04 tetrahedra linked to a continu- 
ous three dimensional network^. Numerical simulations 
in the framework of classical molecular dynamics (MD) 
has allowed in the past to clarify several issues connected 
with the structural properties of this material at the mi- 
croscopic level. Basis of this approach are interaction 
potentials fitted to ab initio and experimental param- 
eters. Density Functional Theory (DFT)^ provides in- 
stead a description of both the electronic and ionic de- 
grees of freedom and represents a powerful tool to study 
disordered systems. Indeed, this approach allows to im- 
prove the description of the structure and gives access to 
the electronic properties of amorphous 5^02 [a-SiO^^- 
Car-Parrinello molecular dynamics (CPMD) simulation 
allows to study dynamical properties of the system within 
DFT— . During the nuclear dynamics, the wavefunction 
remain sufficiently close to the Born-Oppenheimer sur- 
face allowing a faithful description of the electronic struc- 
ture in most standard situation. However, the treatment 
of the dynamics of disordered systems in an ab-inito way 
is extremely costly from a computational point of view. 
This imposes restrictions to the size of the systems that 
can be studied and on the maximum lenght of the trajec- 
tory that can be calculated in a molecular dynamics simu- 
lation. The lenght of the trajectory that can be obtained 
within the Car-Parrinello method (CPMD)^ is of the or- 
der of picoseconds. In order to overcome these limitations 
different simulation strategies have been adopted in the 
past that combine the advantages of classical dynamics 
and the accuracy of an ab initio description. Among 
these, the methodology exemplified in a paper by Gin- 
hoven et al.— has been validated in a very accurate way. 
Through the comparison of several combinations of clas- 
sical and density functional schemes, the authors show 



that small systems of a-Si02 (under 100 atoms) exibhit 
local structural characteristics that are similar to those 
of larger systems^. Benoit et ali^ use classical molecular 
dynamics for annealing small samples at very high tem- 
peratures (above T — SOOOiir), quench the liquid up to 
a glassy transition, and then refine their analysis with a 
CPMD dynamics run followed by CPMD quenching to 
find a minimum energy structure. However, the work 
of Benoit et al. still leaves space for improvement. In 
this paper we propose an enhanced protocol based on 
the same philosophy, but with a few important differ- 
ences. First of all we use an improved classical potential, 
suitable for high temperature phases. Subsequently, our 
CPMD quenching protocol, particularly suited for the 
amorphous samples of small size shows to be more effi- 
cient than the one used in the past. Finally, we provide an 
interesting insight into the ring analysis by doing energet- 
ical considerations that should help to assess the validity 
of amorphous samples generated with this method. This 
work shows how to efficently generate small systems that 
still well represent the properties of infinite sized sam- 
ples. The procedure described in this paper reduces the 
amount of CPU time required for the generation of a- 
Si02 by about a factor of five with respect tc^ without 
compromising the final results. The properties of these 
system models are in good agreement with fully ab initio 
model systems and with available experimental data. 



II. 



COMPUTATIONAL SETUP 



Our starting model was a slightly deformed cell of 72 
atoms in the solid a-quartz phase. We adopted periodic 
boundary conditions (PBC) in all directions in order to 
mimic the bulk behavior. The PBC scheme allows to re- 
duce finite size effects as well as surface effects and to 
get realistic properties with a small simulation cell. The 
model of the amorphous phase a-Si02 was obtained by 
first melting the quartz cells and then quenching from the 
melt using classical molecular dynamics with step cooling 
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TABLE I. Parameters of BKS potential. 

A,,(eV) C,j(eVi^) charges 



O-O 
Si-0 
Si-Si 



1388.7730 
18003.7572 



2.76000 
4.87318 



175.0000 
133.5381 



-1.2 
2.4 



TABLE IL Correction parameters to BKS potential. 



charges 



O-O 
Si-O 
Si-Si 



180.00 
20.00 



24.00 
6.00 



is. 



= -1.2 
= 2.4 



from T = 8500 to 300 K at a quenching rate of 1.6 x 10^^ 
K/s (see section III). We then relaxed the obtained con- 
figurations using the formalism of DFT. 
The classical MD calculations were performed using the 
MOLDY package^. The long range electrostatic forces 
were treated using the Ewald technique^. We classically 
modeled the interaction between atoms using the em- 
pirical potential introduced by van Beest et al. (BKS)^ 
which has been shown to properly describe the structural 
and dynamical properties of amorphous silica. This po- 
tential is a two body interaction consisting of a Coulomb 
term and a short range term cast in the usual Buck- 
ingham form. Only two short range interactions are 
taken into account: the Si-0 interaction which describes 
the silica bond and the 0-0 non bonded interaction 
which causes the tetrahedal arrangement of oxygen atoms 
around the silicon atom. At short distances the BKS po- 
tential diverges attractively limiting its use in molecular 
dynamics runs at high temperature (in our case above 
3700 K). In order to avoid this divergence at short dis- 
tances we add a repulsive terniiS,. This modification is 
necessary to describe the high temperature liquid where 
ions of opposite charges would otherwise approach each 
other very closely and be trapped in the well of the 
Coulomb potential. The functional form of the poten- 
tial adopted in the classical simulations is therefore: 



Atteq r.i 



D 



Ei. 



(1) 



where the first three terms represent the standard BKS 
potential and the last two terms the correction intro- 
duced to avoid the divergence at short distances. The 
force field parameters for the BKS potential are reported 
in Table 1. In Table 11 we report the coefficients of the cor- 
rection terms applied and in Fig. [l]we plot the original 
and the modified BKS potential. These additional terms 
do not change the physics of crystalline and amorphous 
phases. 



DFT calculations were performed using the CPMDH 
software. The electronic structure has been treated using 
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FIG. 1. Plot of the interatomic potential used in the classical 
MD simulation. In the upper graph we show the original BKS 
potential describing the interaction between O-O atoms and 
the modified version implemented in the simulations. In the 
middle graph we plot the BKS and the modified potential 
describing the Si-O interaction. In the lower graph we report 
the plot of the potential for the Si-Si interaction. In this case 
we have only Coulomb interactions and no short range terms. 



plane waves (PW) basis sets and pseudopotentials. The 
use of PW has several advantages. PW call for the use 
of PBC and the atomic forces can be computed via the 
Hellman-Feynman theorem without evaluating the Pulay 
contributionsi^. In addition the convergence of total en- 
ergy and force calculations can be controlled by a single 
parameter (kinetic energy cutoff) and improved to arbi- 
trary accuracy. We made use of the Generalized Gra- 
dient approximation (GGA)^"^ to describe the exchange- 
correlation energy. Core-valence interactions have been 
described through norm-conserving pseudopotentials for 
both oxygen and silicon atoms^^. We used a plane wave 
basis set defined by an energy cutoff of 70 Ry for plane 
waves and 280 Ry for the electronic density which has 
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been shown to be large enough to insure the convergence 
of the energy, in some prehminary tests on a-quartz and 
/3-cristobahte. The Brihouin zone has been sampled at 
the r point only^^. All the Density Functional calcula- 
tions have been performed with PBC in all three spa- 
tial directions and keeping the volume of the supercell 
costant. 
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Time [ps] 



III. MODEL SYSTEMS AND QUENCHING 
PROTOCOL PROCEDURE 



FIG. 2. Quenching protocol used in this work for the gener- 
ation of the amorphous Si02 model systems. 



An orthorombic simulation cell of crystalline a-quartz 
containing 72 atoms of Oxygen and Slicon (24 Si02 units) 
has been rescaled in the x, y and z direction making it 
cubic and matching the amorphous Si02 experimental 
density of 2.20 gr/cm-^. The so obtained configuration 
has been used as the starting supercell for the simula- 
tions. The size of this periodically repeated cubic super- 
cell was 10.29 A. 

In the classical MD simulations the equations of motion 
for the ions were integrated using a modified version of 
the Beeman algorithm as implemented in MOLDY^ with 
a timestep of 1.6 fs. In total seven amorphous silica 
model systems were generated. During the simulations 
all the atoms within the supercell have been allowed to 
move while the volume of the supercell has been kept 
fixed. The starting configurations have been heated up 
rapidly to T = 8500 K leading to the formation of liq- 
uid silica. We subsequently performed a 25 ps costant 
temperature run at T = 8500 K in order to equilibrate 
these high temperature liquids. We then quenched the 
model systems to T = 300 K using constant tempera- 
ture MD runs performed at different temperatures, with 
a final averaged quenching rate of 1.6 x 10^^ K/s. The 
quenching protocol is reported in Fig. [2] It is impor- 
tant to notice that the quenching rate should be not too 
high. Very high quenching rates hinder an adequate re- 
laxation, and can eventually freeze the concentrations of 
small rings (see section IV) to values much higher than 
those experimentally expected. During all the costant 
temperature runs the velocities of the particles were pe- 
riodically rescaled assuming a Maxwell-Boltzmann distri- 
bution at the target temperature. We noticed from the 
mean square diplacements of the Si and O atoms, that 
the freezing temperature of these model systems can be 
located around T = 3300 K, which is much larger than 
the experimental data. The resulting structures were 
perfectly ordered with every silicon atom fourfold and 
every oxygen atom twofold coordinated. The final con- 
figurations were equilibrated classically at T = 300 K 
for about one nanosecond. We finally performed a 0.5 
picosecond microcanonical Car-Parrinello molecular dy- 
namics run on the equilibrated systems (T = 300 K) 
in order to get configurations for structural analysis and 
comparison with the classical model systems. The anal- 
ysis of the structural properties is reported in section VI. 



IV. AB INITIO STRUCTURAL RELAXATION 

The configurations obtained as described in the pre- 
vious section have been further optimized by explicitly 
treating the electronic degrees of freedom with DFT. 
DFT combined with Molecular Dynamics offers several 
methods for structure relaxation and global optimization. 
In this section we describe the effects of geometry opti- 
mization (G.O.) and Car-Parrinello simulated annealing 
(CP annealing) of the model systems. G.O. and CP an- 
nealing represent two different ways to reach the equi- 
librium configuration of a given system. In a G.O. one 
starts with a given structure characterized by the nuclear 
coordinates. By evaluating the forces for that set of nu- 
clear coordinates one gains information on how to shift 
the atomic positions in order to reduce the total energy. 
The procedure is repeated until the structure does not 
change any further within a predefined tolerance. An 
alternative approach is the CP annealing. Within this 
method the electronic and nuclear degrees of freedom 
are optimized simultaneously, starting the calculation at 
a given temperature. Then the temperature is reduced 
gradually and the structure should get trapped in that 
of the global energy minimum with a higher probability 
and not in a local one^^ . This method has been shown to 
be successful for large systems where the energy surface 
is extremely frustrated. 

When switching from the Classical MD treatment to the 
ab-initio description, a common way to treat the system 
is to first perfom a sufficentcly long constant temperature 
(NVT) run, say at T = 300 K for equilibration, and then 
anneal to a T = K configuration to reach the ground 
state energy of the systeni^iii. This procedure costs a lot 
of CPU time (ii^iH larger systems have been studied). 
We noticed that, at least for these small systems, a G.O. 
and a CP annealing lead to the same final configuration, 
but the CPU time required for the G.O. to reach conver- 
gence is one fifth of the time needed by a CP annealing 
(see Fig. [3]). In addition during the optimization or the 
annealing no substantial change occurs on the topology 
acquired at the end of the classical quenching protocol. 
The root mean square displacement (RMSD) between the 
final configurations obtained within G.O. and CP anneal- 
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FIG. 3. Total energy of the system versus the CPU time in 
seconds with four processors. The upper graph shows the 
behaviour of the system in a geometry optimization, and the 
lower graph shows the behaviour in a Car-Parrinello simulated 
annealing. 



TABLE III. CPU time in seconds with four processors needed 
to relax the system to its ground state as a function of the 
number of time steps of the initial CP run. 



Step 





100 


200 


300 


500 


700 


2000 


CP (s) 





2314 


4628 


6943 


11572 


16201 


46289 


G.O. (s) 


9731 


9340 


10059 


9602 


6819 


5779 





Total (s) 


9731 


11654 


14687 


16545 


18391 


21980 


46289 



ing methods has been evaluated to be on average 0.012 A 
while the mean difference in energy is 0.002 eV. The same 
convergence criteria (10~^a.u. for the energy and 10~^ 
a.u.-A^-'^ on the gradients) have been adopted in both 
procedures. In the G.O. we used the limited-memory 
gpQgis method as implemented in CPMD. In the CP 
annealing the parameters that can be controlled are the 
fictitious or inertia mass /x, the scaling factors for ions and 
electrons respectevely and 7e and the time stepi^r^i. 
We choose /x=700 a.u., 7Ar=0.99, 7e=0.99 and a time 
step of 0.12 fs. We further investigated the relaxation 
process by applying a combination of CP annealing and 
G.O. to the model systems. In table IIIII we report the 
CPU time needed to reach the convergence by applyng 
first a CP annealing for a certain number of time steps 
followed by a G.O. As can be inferred from the table the 
CPU time is minimum when no annealing is performed 
and the structure is simply geometry optimized. 




10000 20000 30000 40000 50000 
CPU time 

FIG. 4. Temperature behavour of the system during the same 
Car-Parrinello simulated annealing represented in the lower 
graph of Fig. |3l 



CRITERIA FOR SELECTING THE MODEL 
SYSTEMS AND RING ANALYSIS 



One way to characterize Si02 networks is by studying 
the statistics of closed rings of bonded atoma^^. An N 
ring is described as a closed Si-O-Si-0- ■ ■ chain with N 
silicon atoms^"^. Small rings {N=2 or 3) are geometrically 
strained structures and result in an energy penalty^*. 
During the cooling protocol unrealistic geometry struc- 
tures can occur. Two reasons for the appearence of such 
unphysical structures are the limited size of the computa- 
tional box (72 atoms) and the applied PBC. In addition 
a too high quenching rate can prevent an appropriate 
relaxation of the model systems and freezes the concen- 
tration of small rings^. Unrealistic model systems exhibit 
the presence of large voids and compressed regions within 
the simulation box or concentrations of small membered- 
rings much higher than the experimental ones^. While 
large voids can be consistent with the structure of amor- 
phous silica, the high density of compressed regions of 
the cell is not. For these reasons we generated several 
model systems of a-Si02 and selected the ones which 
satisfied all the criteria of a realistic structure. A pri- 
mary criterion for selecting the systems was to compute 
the total energy and eliminate the model systems with 
higher energies corresponding to higher stresses. As a 
consequence of this criterion systems which exibihited a 
too high concentration of small rings, which are strained 
structures, have been eliminated. Among seven gener- 
ated systems only one was eliminated because it showed 
a too high concentration of small rings confirming that 
the quenching protocol adopted was slow enough for an 
appropiate relaxation of the amorphous configurations. 
The three-membered rings are quasiplanar as can be de- 
duced from the sum over all bond angles in the rings that 
average to 702°, while the ideal value is 720°. The aver- 
age Si-O-Si angle in these rings is 131.80° which is smaller 
than the average of 147° of the model systems. The av- 
erage Si-0 bond lenghts is 1.642 A in agreement witb2i. 
Four-membered rings do not exhibit preferences for pla- 
narity; the sum over all bond angles in the rings give an 
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TABLE IV. In the second column AE is the dilfercncc of the 
total energy of the model systems to the one with the lowest 
energy. In the other columns number of n-rings in the model 
systems. 



models AE(eV) 2-R 3-R j-R 5-R 6-R 7-R 

model- 1 ri5 1 1 1 1 T" 

model-2 1.65 3 1 2 1 

model-3 0.7 1 1 3 2 

model-4 1 1 2 1 

model-5 0.44 2 1 2 1 

model-6 1 2 1 3 2 1 

model-7 1.2 12 14 3 1 



3-R 

1 
1 

1 





angle of 976.88° (ideal 1080°) while the average Si-O-Si 
angle is 137.02°, which is smaller than the reference an- 
gles inferred from NMR measurements (145° — 150°'^ ) or 
from x-ray diffraction experiments (144°^^). 
In table IIVI we report the differences AE of the total 
energy of the model systems to the one with the lowest 
energy, which we called model-4. In the other columns we 
give a description of the systems in terms of ring statis- 
tics. As can be seen from the table the systems with 
higher energies and stress are that ones with a higher 
concentration of small rings, showing that the presence 
of small rings is the main cause of stress in the system. 
We estimated an upper bound fo AE of 1 eV per small 
rings (2,3 membered rings) which is in agreement with—, 
considering the fact that in additon to the penalty en- 
ergies coming from the strained rings other sources for 
stress can be the fact that during all the simulations the 
volume of the computational box has been kept fixed. 



the angle distibutions is a small shift of the mean value of 
the intra-tetrahedral Si-O-Si angle. We estimated for this 
angle a mean value of 152° ±11° for the classical aver- 
aged configurations and 146° ± 6° for the Car-Parrinello 
microcanonical run (NVE) (the experimental value being 
140°-150°^'^). We noticed that this change in angle takes 
place during the first steps of the Car-Parrinello run cor- 
responding to an extremely short sampling time (« 0.07 
ps). 

We have then calculated the averaged pair correlation 
functions gap{r), where a, (3— Si, O of the systems at 
T = 300 K for the Classical and Car-Parrinello runs 
(Fig. |6]). Within the classical treatment we obtained 
an averaged distance of I.6IA for Si-0, 3.12A for Si-Si, 
and 2.66A for 0-0, in very good agreement with the ex- 
perimental data^^. For the Car-Parrinello treatment the 
averaged bond lenght Si-0 is found to be 1.65 A, which is 
slightly (3%) larger than the experimental value (I.6IA), 
in agreement with a general tendency of the GGA^'*. Us- 
ing the same computational setup we found a similar 
overstimation in a preliminary test made on a-quartz. 
For the other peaks we found 3.18A for Si-Si and 2.68A 
for the 0-0 bond lenght. 

From this analysis we concluded that our classical sys- 
tems were well equilibrated. From the pair correla- 
tion functions we computed the static structure factor 
for the Car-Parrinello run which can be compared with 
experiments^^. Experimentally the static structure fac- 
tor can be obtained from neutron diffraction. In the sim- 
ulations one possibility is to compute the static structure 
factor S{q) by its relation to the pair correlation function 



VI. STRUCTURAL PROPERTIES 
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Among all the generated Si02 model systems we se- 
lected the model-4 which best satisfied our criteria for a 
detailed study. This choice has been performed on the 
basis of the ring analysis of the previous section, but 
further analysis has to be performed to validate this fi- 
nite size sample. The most natural way to further esti- 
mate the quality of the generated amorphous model sys- 
tems is to analyze their structural characteristics such as 
bond lenghts and angle distributions, radial distribution 
functions, static structure factors and to compare them 
to available experimental data. We evaluated all these 
properties for the system called model-4. At T = 300 
K we have calculated the time-averaged distributions of 
the intra-tetrahedral 0-Si-O and inter-tetrahedral Si-O- 
Si angles for the classical and quantum MD simulations 
using the same starting configuration (Fig. [5]). The intra- 
tetrahedral angles 0-Si-O stay close to the experimental 
value of 109.4° for both simulations: 109.3° ± 5 for the 
classical and 109.5° ± 7 for the CPMD. The only dif- 
ference between the classical and quantum simulation of 



(3) 



The integral in Eq.([21) has been evaluated by using 
the Filon's method'''^. In equation ([2]) p stands for the 
density. In Eq. ^ the gap{r) are the pair correlation 
functions (ck,/3 — Si,0) while Ca and cp stand for the 
concentrations of the two species and ba and bp are their 
scattering lenghts. We chose and 6/3 to be 4.149 fm 
for Si and 5.803 fm for O^^. The integration from to 
00 is performed as a summation from zero to half of 
the box size. The first sharp diffraction peak (FSDP) 
corresponds to 1.53A~^ which is in good agreement with 
the eperimental data (1.52A~^)^. The other peaks are 
at 2.8, 5.19, and 7.83A~^, respectevely. 
We finally studied the electronic density of states of 
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FIG. 5. Averaged angle distribution for O-Si-O (upper graph 
and Si-O-Si (lower graph). Solid lines refer to Car-Parrinello 
MD and dashed lines refer to Classical MD simulations. 



the system. We evaluated the Kohn-Sham density 
of states (Fig. [S]). The density is in agreement with 
other DFT calculations^!^. The calculated band gap 
is 5.7 eV and it is understimated (experimental ^ 
9 eV) as usual in DFT due to the discontinuity of 
derivatives of the exchange-correlation functional^. 
The DOS of amorphous and crystalline Si02 are 
very similar, because the short range order is similar 
for both structures. For this reason we can identify 
the states at about -20 eV as oxygen 2s states, the 
states from -10 eV to -4 eV as bonding states between 
silicon sp^ hybrids and oxygen 2p orbitals and above 
-4 eV as oxygen 2p nonbonding orbitals. The lowest 
conduction band are states with antibonding character—. 
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FIG. 6. Simulated pair correlation functions for Si-O (upper 
graph), Si-Si (middle graph) and 0-0 (lower graph). The pair 
correlation functions are obtained from NVE simulations at 
300 K. The solid lines refer to Car-Parrinello runs while the 
dashed lines refer to Classical MD runs. 



Static neutron structure factor S(q) 




VII. 



CONCLUSIONS 



Using classical molecular dynamics and the BKS em- 
pirical potential we generated seven model systems of 72 
atoms of a-Si02 at 300K. We further optimized the clas- 
sical equilibrated configurations with plane wave based 
density functional theory and GGA approximation to the 
exchange correlation energy. In particular within ab ini- 
tio treatment we showed that both geometry optimiza- 
tion and Car-Parrinello annealing applied to our systems 
lead to the same final configurations but the CPU time 
required for the geometry optimization to reach conver- 
gence is one fifth of the time needed by a Car-Parrinello 
annealing. The topology acquired by the amorphous sil- 
ica model systems at the end of the classical quench- 
ing protocol did not change when treated in an ab initio 
way. The structural properties of the amorphous glass, 
like angle distributions, pair correlation functions, static 
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FIG. 7. S(q) at K computed with CPMD (upper graph), 
and experimental results from neutron diffraction of normal 
and densifled vitreous Si02 (from Ref.-^) (lower graph). 
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FIG. 8. Electronic density of states of amorphous Si02 at OK. 



structure factor and electronic density of states have been 
accurately studied. The properties of rings in model sys- 
tems have been accurately analyzed by energetical con- 



siderations. We found very similar results for both the 
description (classical and Car-ParrincUo) and in good 
agreement with experimental data. 
The generated model systems will be used to study de- 
fects, and in particular the mechanism for bond weaken- 
ing and electron localization in 5^02. 
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